One of the most common equations that start your adventure with the finite element method is the heat transport equation. Our unknown function \( u \) here is a function that assigns a temperature value to various points in the area in which we are solving our problem. In other words, \( u(x,y) \) represents the temperature value at point \( (x,y) \) of the domain \( \Omega \).
The temperature distribution function sought is defined as follows: \( \Omega \ni (x,y) \rightarrow u(x,y)\in R \).
The strong wording for the heat transport problem is as follows \( -\Delta u=f \)
In this equation \( -\Delta u \) stands for the so-called Laplasjan, i.e. the sum of the second-degree partial derivatives, and to simplify the notation, we omitted the variables \( (x,y) \).
\( \Delta u (x,y) = \frac{\partial^2 u(x,y) } { \partial x^2} + \frac{\partial^2 u(x,y) }{\partial y^2 } \)
while \( f(x,y) \) means the function representing the heat source. In those parts of the domain where the function \( f(x,y) \) is positive, there is "heating", while in those parts of the domain where the function \( f(x,y) \) is negative, there is "cooling".

The weak wording is obtained as follows. We integrate and multiply our equation by selected functions \( v(x,y) \) called test functions which we will use to average our equation in the area where these functions are defined.

In a similar way to the bitmap projection problem, any choice of test function \( v(x,y) \) for averaging our problem in the area where the test function is defined gives us one equation. Various test function selections \( v(x,y) \) gives different equations ( 1 ).
If we write our formula for partial derivatives, we get two separate integrals \( \int_{\Omega} \Delta u(x,y)v(x,y) dxdy= \int_{\Omega} \left( \frac{\partial^2 u(x,y) }{\partial x^2 } + \frac{\partial^2 u(x,y) }{\partial y^2 } \right) v(x,y) dxdy \)
The problem that we want to solve looks like this: We are looking for the temperature distribution, a function \( \Omega \ni (x,y) \leftarrow u(x,y) \in R \) such that
\( \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial x^2} v(x,y) dxdy + \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial y^2} v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)
for any test functions \( \Omega \ni (x,y) \rightarrow v(x,y) \in R \) (of course our temperature distribution and the testing function must be such that it is possible to calculate the integrals, as we know, not all functions can be used to calculate integrals, some functions give integrals that can be infinite),
\( a(u,v)=l(v) \\ a(u,v)=\int_{\Omega } \frac{\partial^2 u(x,y) }{\partial x^2 } v(x,y) dxdy + \int_{\Omega } \frac{\partial^2 u(x,y) }{\partial y^2 } v(x,y) dxdy \\ l(v)=\int_{\Omega} f(x,y) v(x,y) dxdy \)
Now note that if we were to now apply our B-spline functions introduced in the previous section to the temperature field approximation, \( u(x,y) = \sum_{i=1,j=1}^{N_x,N_y} u_{i,j} B^x_{i}(x) B^y_{j}(y) \)
then we would have to be able to compute the second order derivatives of the B spline function.
In order for the second-degree derivatives to be computable, the basis functions used to approximate the temperature field must be smooth. B-splines are great for calculating derivatives of these. However, in order for these derivatives not to be zero (so that our system of equations, which is built similarly to bitmap approximation, does not end up equal to zero), the B-splines used here must be higher-order functions (at least of the class \( C^2 \), i.e. they must be at least 3rd degree B-splines).
If we wish to use higher order polynomials, we can now insert the integrals of our equation into the problem matrix, i.e. replace the bitmap projection problem matrix with our matrix.
Now, in place of functions \( u \) we insert an approximation using the B-spline function \( u = \sum_{i=1,j=1}^{N_x,N_y} u_{i,j} B^x_{i} B^y_{j} \), and in place of the test functions we also put B-spline function \( B^x_{k} B^y_{l}, k=1,...,N_x; j=1,...,N_y \). Similarly to the bitmap projection, we get a system of equations
\( \begin{bmatrix} a({B^x_{1,p}B^y_{1,p}},B^x_{1,p}B^y_{1,p}) & a({B^x_{1,p}B^y_{1,p}},{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{1,p}B^y_{1,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ a({B^x_{2,p}B^y_{1,p}},{B^x_{1,p}B^y_{1,p}}) & a({B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{2,p}B^y_{1,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ \vdots & \vdots & \vdots & \vdots \\ a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{1,p}B^y_{1,p}}) & a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y}\\ \end{bmatrix} \\ = \begin{bmatrix} l( B^x_1(x)*B^y_1(y)) \\ l(B^x_1(x)*B^y_2(y)) \\ \vdots \\ l(B^x_{N_x}(x)*B^y_{N_y}(y)) \\ \end{bmatrix} \)
where after taking into account the formulas on \( a \) i \( l \) we have

Of course, we must remember that
\( \int_{\Omega} \Delta u(x,y)v(x,y) dxdy= \int_{\Omega} \left( \frac{\partial^2 u(x,y) }{ \partial x^2 } + \frac{\partial^2 u(x,y) }{ \partial y^2 } \right) v(x,y) dxdy \)
Now, if we generated the appropriate integrals and run one of the matrix factorization algorithms described in Chapter 5, it would turn out that the matrix is singular and the system of equations cannot be solved! Why would that happen? Because when we solve the heat transport problem and define the equation satisfied in the area \( \Omega \) we need to specify the so-called boundary conditions, i.e. describe how the temperature behaves on the boundary of the area.
Whenever we want to formulate our engineering problem with partial differential equations, for example the problem of heat transport, over some area
\( \Omega \), we must first of all give the partial differential equation that describes the modeled phenomenon inside the region \( \Omega \) Additionally, we must describe the so-called boundary conditions, i.e. provide some additional equations that describe the behavior of our phenomenon at the boundary of the area.
Let us imagine that we want to estimate the temperature inside the room. Of course, we need to specify what partial differential equation describes the physical phenomenon of heat flow inside the room. This is of course the already known equation (strong formulation) with Laplasian
\( \Delta u(x,y)=f(x,y) \). If we give the equation itself only for the center of the room, but we do not take into account what is happening on the walls of the room, we have potentially an infinite number of possible solutions. So, our system of equations cannot be solved without specifying what is happening on the boundary. The correct definition of the boundary conditions usually gives us one possible solution.
Our area \( \Omega \) can represent the room interior. For the sake of simplicity, let's assume that we are working in two dimensions (that our room is flat and two-dimensional). Let us assume that this room is L-shaped and its inner walls are adjacent to a radiator on the boundary of which the temperature is zero, and the remaining walls have windows through which external heat reaches us. Then the temperature on the inner walls of the room will be zero. We introduce the so-called Dirichlet boundary condition on the part of the boundary of the area (on the internal walls on which we measured the temperature). For these walls (part of the boundary \( \partial \Omega \)) we denote symbols \( \Gamma_D \) and we give the temperature values on this part of the boundary (on the walls of the area)
\( u(x,y) = T(x,y) \textrm{ dla }(x,y)\in \Gamma_D \)
where \( T(x,y) \) this is the temperature measured with a thermometer or a thermal imager. This situation is illustrated in Fig. 1.

There may also be windows on a certain part of the boundary of the room, and outside the windows it can be hot (if it is summer) or cold (if it is winter) and the temperature around the glass can rise (in summer) or decrease (in winter) at any time. Then, by using a thermometer and a watch (or a ‘magic’ thermal imaging camera), we can measure the speed with which the temperature changes next to the glass.
In this case, we introduce the so-called Neumann boundary condition on the part of the boundary of the area (on the windows where we measured the speed of temperature changes). This part of the boundary (i.e., the window) is marked with the symbol \( \Gamma_N \) and we determine the speed of temperature changes
\( \frac{\partial u(x,y)}{\partial n} = g(x,y) \textrm{ dla }(x,y)\in \Gamma_N \)
where \( g(x,y) \) is the heat flux on the boundary of the area.
What does it mean here \( \frac{\partial u }{\partial n } \)? This indicates the speed of temperature change towards the boundary. From a mathematical point of view \( \frac{\partial u }{\partial n }=\frac{\partial u }{\partial x }n_1+\frac{\partial u }{\partial y }n_2 \) where \( n=(n_1,n_2) \) is the vector perpendicular to the boundary of the area \( \Omega \), directed outwards, i.e. if the boundary is a wall perpendicular to the axis \( x \) to \( n=(1,0) \), if the boundary is the bottom edge of the room, is \( n=(0,-1) \) itd. Sample vectors \( n=(n_1,n_2) \) shown in Fig. 1.
How would we modify our equation system ( 2 ) to include zero temperature on part of the boundary \( \Gamma_D \) and the speed of heat flow through the part of boundary \( \Gamma_N \).
For this purpose, we need to specify what our element group and our B-spline basis functions will look like spanned over the elements, giving a vector of nodes along the axis \( x \) and the knot vector along \( y \). TThe simplest solution is to divide our area into four square sub-areas (four groups of elements) and glue them together (repeating the nodes on the glue interface so that the basis functions at this point have C0 continuity). Here we refer to the third chapter, which describes the way of defining basis functions with the help of node vectors and B-spline polynomials.
Along the axis \( x \) we introduce a vector of nodes [0 0 0 1 2 2 3 4 4 4], similarly along the y axis we introduce a vector of nodes [0 0 0 1 2 2 3 4 4 4].
We repeat the middle node to reduce the continuity of the B-spline in the middle of the area. This will allow us to correctly enforce the Dirichlet boundary condition on one area quadrant. Note that the B-spline functions are smooth and span many elements, so it is impossible to force the values of the B-spline function at some point without reducing continuity (without repeating nodes in node vectors).
The resulting basis functions are shown in Fig. 2. We obtained two bases of one-dimensional B-spline functions \( \{ B_{i,2}(x) \}_{i=1,...,5 } \) and \( \{B_{j,2}(y)\}_{j=1,...,5 } \). Then we will create tensor products from them \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,5 \).

Note that our area \( \Omega \) is stretched over a square \( [-1,1]^2 \). But our 5 * 5 = 25 basis functions are spread over a square \( [0,2]^2 \) resulting from the definition of the node vector [0 0 1 2 2 3 4 4].
So we need to "remap" our physical area \( [-1,1]^2 \) to a group of elements on which our functions are defined B-spline
\( [-1,1]^2 \ni (x,y) \rightarrow \psi(x,y)= \left(2(x+1),2(y+1)\right) \)
We extend our functions to the area \( [-1,1]^2 \) assuming that outside our region they are equal to zero, and we change the variables in the block integral \( [-1,1]^2 \) per block \( [0,4]^2 \)
Then we remap all the integrals to our area where we span the integrals \( a(u,v)=\int_{\Omega} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) dxdy + \int_{\Omega} \frac{\partial^2 u(x,y) }{ \partial y^2 } v(x,y) dxdy= \int_{[-1,1]^2} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) dxdy + \int_{[-1,1]^2 } \frac{\partial^2 u(x,y) }{ \partial y^2 } v(x,y) dxdy= \\ \int_{[0,4]^2} \frac{\partial^2 u(x,y) }{\partial x^2 } v(x,y) |Jac(\psi^{-1})| dxdy + \int_{[0,4]^2} \frac{\partial^2 u(x,y)}{ \partial y^2 } v(x,y) |Jac(\psi^{-1})| dxdy \\ l(v)=\int_{\Omega} f(x,y) v(x,y) dxdy = \\= \int_{[-1,1]^2} f(x,y) v(x,y) dxdy = \int_{[0,4]^2} f(x,y) v(x,y) |Jac(\psi^{-1})| dxdy \)
In the above formulas for the integrals, we put the absolute value of the Jacobian transformation of the change of variables from our area Ω=[−1,1]2 on the area \( \Omega=[-1,1]^2 \) on the area \( [0,4]^2 \) on which we unbuttoned the knot vectors. It is calculated like this. We know the mapping \( \psi:\Omega \rightarrow [0,4]^2 \) dane \( \psi(x,y)=\left(2(x+1),2(y+1)\right) \). This means that new variables in the area \( [0,4]^2 \) (which we will mark for a moment \( (\xi,\eta) \) data are formulas \( (\xi,\eta)= \left(2(x+1),2(y+1)\right) \). If I would like to calculate the variables \( (x,y) \) in the area \( [-1,1]^2 \) I just have to convert it
\( 2(x+1)=\xi \) so \( x=\xi/2-1 \), and \( 2(y+1)=\eta \), so \( y=\eta/2-1 \).
So \( \psi^{-1}(\xi,\eta)=(\xi/2-1,\eta/2-1) \). For convenience, let's change the variables back to symbols \( (x,y) \) czyli \( \psi^{-1}(\xi,\eta)=(x/2-1,y/2-1) \). Now, we need to count the partial derivatives \( \frac{\partial(x/2-1)}{\partial x }=1/2 \), \( \frac{\partial(x/2-1) }{\partial y }=0 \),
\( \frac{\partial(y/2-1) }{\partial x }=0 \),
and \( \frac{\partial(y/2-1) }{\partial y }=1/2 \). Our Jacobian integral is the determinant of the matrix of partial derivatives
\( jac(\psi^{-1})=|\begin{bmatrix} 1/2 & 0 \\ 0 & 1/2 \end{bmatrix}| = |1/2*1/2|=|1/4|=1/4 \)
that is, our system of equations
\( a(u,v)=\int_{[0,4]^2} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) 1/4 dxdy + \int_{ [0,4]^2 } \frac{ \partial^2 u(x,y) }{\partial y^2 } v(x,y) 1/4 dxdy \\ l(v)=\int_{[0,4]^2 } f(x,y) v(x,y) 1/4 dxdy \)
Now, for the B-spline functions spanning on the edge of the Dirichlet area, we can force the Dirichlet condition in such a way that we zero the appropriate lines and right sides of the system of equations, and enter 1s on the diagonal.
In other words, for the relevant lines corresponding to the B-splines \( i,j \) which are on the Dirichlet boundary
\( \int{\Delta(B^x_{i,p }B^y_{j,p } ) }{B^x_{1,p } B^y_{1,p } }u_{1,1 } + \cdots + \int{ \Delta(B^x_{i,p }B^y_{j,p}) }{ B^x_{2,p }B^y_{1,p} }u_{i,j }+ \cdots +\int{ \Delta( B^x_{i,p }B^y_{j,p }) }{B^x_{N_x,p }B^y_{N_y,p } }u_{N_x,N_y } = \\= \int f(x,y)B^x_{ i,p }B^y_{j,p } dxdy \\ \)
we replace them with \( 1*u_{i,j }=0 \)
But how to express the Neuman condition? Unfortunately, this is not easily done without integration by parts.
So we have to apply the formula for integration by parts. Since our problem is two-dimensional, the formula is relatively complicated. We will not introduce it here because it would require taking a course in mathematical analysis. \( \int_{\Omega} \Delta u (x,y) v(x,y) dxdy =- \int_{\Omega} \nabla u (x,y) \cdot \nabla v(x,y) dxdy + \int_{\partial \Omega} \frac{\partial u }{\partial n } v dS \)
where \( \nabla u(x,y) = (\frac{\partial u}{\partial x }, \frac{\partial u }{\partial y }) \) denotes gradient \( u \) that is a vector of partial derivatives of \( u \), similarly \( \nabla v(x,y) = (\frac{\partial v}{\partial x }, \frac{\partial v }{\partial y }) \)
stands for gradient \( v \) that is, a vector of partial derivatives of a function \( v \), while
\( \cdot \) denotes the dot product of their vectors, i.e.
\( (\frac{\partial u }{\partial x }, \frac{\partial u }{\partial y }) \cdot (\frac{\partial v }{\partial x }, \frac{\partial v }{\partial y }) = \frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \).
So our entire equation \( -\int_{\Omega} \frac{\partial^2 u(x,y)}{\partial x^2 } v(x,y) dxdy - \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial y^2} v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)
after integration by parts it turns into \( \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \right) dxdy= \int_{\partial \Omega } \frac{\partial u }{\partial n } v dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)
What does it mean here \( \frac{\partial u }{\partial n } \)? It indicates the speed of temperature change towards the boundary. From a mathematical point of view \( \frac{\partial u }{\partial n}=\frac{\partial u }{\partial x }n_1+\frac{\partial u }{\partial y }n_2 \) where \( n=(n_1,n_2) \) is the vector perpendicular to the edge of the area \( \Omega \), that is, for example, if the edge is a vertical face perpendicular to the axis \( x \) to \( n=(1,0) \), if the edge is the floor of a room \( n=(0,-1) \) etc.
The field in which we solve ours problem \( \Omega = [-1,1] \times [-1,1] \ [-1,0] \times [-1,0] \) has the shape of an inverted letter L.
On the inner part of the edge of the area we mark \( \Gamma_D=\{ (x,y): x=0, y\in[-1,0] \textrm{ or } x \in [-1,0], y=0 \} \) we assume the zero temperature expressed by the Dirichlet boundary condition
\( u(x,y) = 0 \textrm{ dla }(x,y)\in \Gamma_D \)
On the rest of the edge of the area we mark \( \Gamma_N=\Gamma \ \Gamma_D \) we introduce the so-called Neumann boundary condition.
This part of the edge (i.e. the window) is marked with the symbol
\( \Gamma_N \) and we determine the speed of temperature changes
\( \frac{\partial u(x,y)}{\partial n} = g(x,y) \textrm{ dla }(x,y)\in \Gamma_N \)
where \( g(x,y) \) is the set speed of temperature changes on the shore.
In conclusion, the edge of our area \( \partial \Omega \) was dicided into two fragments \( \partial \Omega = \Omega_D \cup \Omega_N \) and on the first we measured the temperature, on the second we measured the speed of temperature changes.
So, we can break a fragment of our weak formulation with the edge integral into these three fragments of the border
\( \int_{\partial \Omega } \frac{\partial u(x,y) }{\partial n } v dS = \int_{\partial \Omega_D } \frac{ \partial u(x,y) }{\partial n } v(x,y) dS + \\ \int_{\partial \Omega_N } \frac{\partial u(x,y) }{\partial n } v(x,y) dS \)
What to do next with these three integrals spanning over three parts of the border?
The integral with the Neumann boundary condition is simple, because our Neumann boundary condition implies that
\( \frac{\partial u(x,y) }{\partial n }=g(x,y) \), where \( g(x,y) \) is a given rate of temperature change, i.e. \( \int_{\partial \Omega_N } \frac{ \partial u(x,y) }{\partial n } v(x,y) dS=\int_{\partial \Omega_N } g(x,y)v(x,y) dS \). Why do we have to convert the derivative in the normal direction into a given function in this integral \( g \)? Since we want to set this boundary condition, we want to model in our simulation a situation where the heat really affects the interior of our room through the window.
As for the Dirichlet boundary condition, it requires some additional treatment.
If \( T(x,y)=0 \) in the Dirichlet condition, i.e. \( u(x,y) = 0 \textrm{ dla }(x,y)\in \Gamma_D \) then the matter is simple. We know that the temperature is zero on this edge, so we do not have to calculate the average values at this point of the edge, we do not have to choose such testing functions \( v(x,y) \) which are not zero on the Dirichlet edge. In other words, I can assume that the testing functions \( v(x,y)=0 \) dla \( (x,y) \in \Gamma_D \). In other words, then the integral
\( \int_{\partial \Omega_D } \frac{\partial u(x,y) }{\partial n } v(x,y) dS=0 \).
Our right side
\( l(v) = \int_{\partial \Omega_D } \frac{\partial u(x,y) }{\partial n } v(x,y) dS + \int_{\partial \Omega_N } \frac{\partial u(x,y) }{\partial n } v(x,y) dS \\ +\int_{\Omega } f(x,y) v(x,y) dxdy = \int_{\partial \Omega_N } g(x,y) v(x,y) dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)
So our problem looks like this
\( a(u,v)=l(v) \\ a(u,v) = \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \right)v(x,y) dxdy \\ l(v) = \int_{\partial \Omega_N } g(x,y) v(x,y) dS \nonumber +\int_{\Omega } f(x,y) v(x,y) dxdy \)
Note that since we have integrated by parts, then we can allow ourselves to reduce the degree of the B-spline function (now, our formula contains only first-order flares, and we are satisfied with the class functions \( C^1 \)).
Moreover, the higher the degree of B-splines, the more expensive it is to process. Simply generating the formula for a higher order B spline requires more operations, as does calculating the value of such a B spline. Bear in mind that the system of equations generated for higher order B-splines will be denser, will have fewer zeros, and therefore the factorization cost (the cost of solving such a system of equations) will be higher.
The solution obtained with isogeometric finite element method is presented in Fig. 3.